---
title: "SALTY_Bivariate"
author: "James"
date: "2023-04-28"
output:
  html_document:
    df_print: paged
---

```{r Packages, echo = FALSE} 
require(OpenMx)
library(stringr)
require(umx)
require(haven)
require(broom)
require(tidyverse)
require(stargazer)
require(eeptools)
require(ggplot2)
require(eeptools)
set.seed(1234)
```

```{r Pre-Processing}

#Loading data
STR <- data.frame(read_dta("1 - Data/cohort_combined.dta"))
Salty = data.frame(read_dta("1 - Data/slty_map.dta"))

combined = merge(Salty, STR, by.x = "LopNr", by.y = "LopNr")

#Removing missing values and incomplete pairs
combined[combined == 99] <- NA
combined[combined == 999] <- NA

#Trust

combined$oldtrust = ifelse(combined$ATTITYD31 == 9, NA, combined$ATTITYD31-1) 

combined$trust = (3-combined$oldtrust)/3

combined$oldtrust = 3-combined$oldtrust

combined$trust_scaled = as.numeric(scale(combined$trust))

#Behavioral

combined <- combined %>% mutate(
                          contactpol = 
                          ifelse(is.na(ATTITYD38_1) == F, yes = 1, no = 0), 
                          contactpub = 
                          ifelse(is.na(ATTITYD38_2) == F, yes = 1, no = 0), 
                          protest = 
                          ifelse(is.na(ATTITYD38_3) == F, yes = 1, no = 0), 
                          boycott = 
                          ifelse(is.na(ATTITYD38_4) == F, yes = 1, no = 0), 
                          finance = 
                          ifelse(is.na(ATTITYD38_5) == F, yes = 1, no = 0), 
                          petition = 
                          ifelse(is.na(ATTITYD38_6) == F, yes = 1, no = 0), 
                          behavioral = (contactpol + contactpub + protest + boycott + finance + petition), 
                          interpersonal_1 = ATTITYD2, 
                          interpersonal_2 = ATTITYD3, 
                          interpersonal_scaled = as.numeric(scale(interpersonal_1 + interpersonal_2)))

#Turnout
Turnout_1970 = data.frame(read_dta("1 - Data/Lev_valdelt_1970.dta"))
Turnout_1994f = data.frame(read_dta("1 - Data/Lev_valdelt_1994f.dta"))
Turnout_1994rkl = data.frame(read_dta("1 - Data/Lev_valdelt_1994rkl.dta"))
Turnout_2010 = data.frame(read_dta("1 - Data/vd10.dta"))
Turnout_2009 = data.frame(read_dta("1 - Data/vd09.dta"))
Turnout_2018 = data.frame(read_dta("1 - Data/Lev_vd18.dta"))
Turnout_2019_EU = data.frame(read_dta("1 - Data/EU2019.dta"))

Turnout_1994 = merge(Turnout_1994f, Turnout_1994rkl, by.x = "LopNr", by.y = "LopNr", all = TRUE)
Turnout_20th = merge(Turnout_1994, Turnout_1970, by.x = "LopNr", by.y = "LopNr", all = TRUE) %>% 
  rename("Ref_1994" = "f", 
         "Ref_1994_Abroad" = "utlandsrost.x",
         "Abroad_1994" = "utlandsrost.y",
         "Nat_1994" = "r.x", 
         "Reg_1994" = "l.x", 
         "Mun_1994" = "k.x", 
         "Nat_1970" = "r.y", 
         "Reg_1970" = "l.y", 
         "Mun_1970" = "k.y", 
         "Legal_1970" = "omyndig", 
         "Foreign_1970" = "utlmedb")

Turnout_0s = merge(Turnout_2009, Turnout_2010, by.x = "LopNr", by.y = "LopNr", all= TRUE) %>% 
  rename("EU_2009" = "e",  
         "Nat_2010" = "r", 
         "Reg_2010" = "l", 
         "Mun_2010" = "k", 
         "Eligible_2010" = "Rostratt.x")

Turnout_10s = merge(Turnout_2018, Turnout_2019_EU, by.x = "LopNr", by.y = "LopNr", all = TRUE) %>% 
  rename("EU_2019" = "Eurost", 
         "Nat_2018" = "rrost", 
         "Reg_2018" = "lrost", 
         "Mun_2018" = "krost", 
         "Eligible_2018" = "Rostratt.x")

Turnout_most = merge(Turnout_20th, Turnout_0s, by.x = "LopNr", by.y. = "LopNr", all = TRUE)

Turnout_all = merge(Turnout_most, Turnout_10s, by.x = "LopNr", by.y = "LopNr", all = TRUE) %>% 
  select(LopNr, EU_2019, Nat_2018, Reg_2018, Mun_2018, Eligible_2018, Eligible_2010, Mun_2010, Reg_2010, Nat_2010, EU_2009, Mun_1970, Reg_1970, Nat_1970, Mun_1994, Reg_1994, Nat_1994, Ref_1994_Abroad, Ref_1994, Abroad_1994, 
         Legal_1970, Foreign_1970)

combined = merge(Turnout_all, combined, by.x = "LopNr", by.y = "LopNr", all.y = TRUE) 

combined = combined %>% 
  mutate(EU_2009 = dplyr::recode(EU_2009, "1" = 0, 
                                 "2" = 1, 
                                 "3" = 2, 
                                 "4" = 1, 
                                 "5" = 1, 
                                 "6" = 1), 
         EU_2009 = na_if(EU_2009, 2), 
         Ref_1994 = dplyr::recode(Ref_1994, "1" = 0, 
                                 "2" = 1, 
                                 "3" = 2, 
                                 "4" = 1, 
                                 "5" = 1, 
                                 "6" = 1), 
         Ref_1994 = na_if(Ref_1994, 2),
         Nat_1970 = case_when(Nat_1970 == "1" ~ 0, 
                              Nat_1970 == "2" ~ 1, 
                              Nat_1970 == "3" ~ NA, 
                              Nat_1970 == "4" ~ 1, 
                              Nat_1970 == "5" ~ 1, 
                              Nat_1970 == "6" ~ 1, 
                              Legal_1970 == "1"|Foreign_1970 == "1" ~ NA), 
         Reg_1970 = case_when(Reg_1970 == "1" ~ 0, 
                              Reg_1970 == "2" ~ 1, 
                              Reg_1970 == "3" ~ NA, 
                              Reg_1970 == "4" ~ 1, 
                              Reg_1970 == "5" ~ 1, 
                              Reg_1970 == "6" ~ 1, 
                              Legal_1970 == "1"|Foreign_1970 == "1" ~ NA),
         Mun_1970 = case_when(Mun_1970 == "1" ~ 0, 
                              Mun_1970 == "2" ~ 1, 
                              Mun_1970 == "3" ~ NA, 
                              Mun_1970 == "4" ~ 1, 
                              Mun_1970 == "5" ~ 1, 
                              Mun_1970 == "6" ~ 1, 
                              Legal_1970 == "1"|Foreign_1970 == "1" ~ NA), 
         Nat_1994 = case_when(Nat_1994 == "1" ~ 0, 
                              Nat_1994 == "2" ~ 1, 
                              Nat_1994 == "3" ~ NA, 
                              Nat_1994 == "4" ~ 1, 
                              Nat_1994 == "5" ~ 1, 
                              Nat_1994 == "6" ~ 1), 
         Reg_1994 = case_when(Reg_1994 == "1" ~ 0, 
                              Reg_1994 == "2" ~ 1, 
                              Reg_1994 == "3" ~ NA, 
                              Reg_1994 == "4" ~ 1, 
                              Reg_1994 == "5" ~ 1, 
                              Reg_1994 == "6" ~ 1), 
         Mun_1994 = case_when(Mun_1994 == "1" ~ 0, 
                              Mun_1994 == "2" ~ 1, 
                              Mun_1994 == "3" ~ NA, 
                              Mun_1994 == "4" ~ 1, 
                              Mun_1994 == "5" ~ 1, 
                              Mun_1994 == "6" ~ 1), 
         Nat_2010 = case_when(Nat_2010 == "1" ~ 0, 
                              Nat_2010 == "2" ~ 1, 
                              Nat_2010 == "3" ~ NA, 
                              Nat_2010 == "4" ~ 1, 
                              Nat_2010 == "5" ~ 1, 
                              Nat_2010 == "6" ~ 1, 
                              Eligible_2010 == "3" ~ NA), 
         Reg_2010 = case_when(Reg_2010 == "1" ~ 0, 
                              Reg_2010 == "2" ~ 1, 
                              Reg_2010 == "3" ~ NA, 
                              Reg_2010 == "4" ~ 1, 
                              Reg_2010 == "5" ~ 1, 
                              Reg_2010 == "6" ~ 1, 
                              Eligible_2010 == "2" ~ NA), 
         Mun_2010 = case_when(Mun_2010 == "1" ~ 0, 
                              Mun_2010 == "2" ~ 1, 
                              Mun_2010 == "3" ~ NA, 
                              Mun_2010 == "4" ~ 1, 
                              Mun_2010 == "5" ~ 1, 
                              Mun_2010 == "6" ~ 1, 
                              Eligible_2010 == "2" ~ NA), 
         Nat_2018 = case_when(Nat_2018 == "1" ~ 1, 
                              Nat_2018 == "0" ~ 0, 
                              
                              Eligible_2018 == "3" ~ NA), 
         Reg_2018 = case_when(Reg_2018 == "1" ~ 1, 
                              Reg_2018 == "0" ~ 0, 
                              
                              Eligible_2018 == "2" ~ NA), 
         Mun_2018 = case_when(Mun_2018 == "1" ~ 1, 
                              Mun_2018 == "0" ~ 0, 
                              
                              Eligible_2018 == "2" ~ NA)) %>% 
  mutate(National_total = rowSums(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970")], na.rm = TRUE), 
         National_eligible = rowSums(!is.na(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970") ])), 
         Regional_total = rowSums(combined[,c("Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970") ], na.rm = TRUE), 
         Regional_eligible = rowSums(!is.na(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970") ])), 
         Municipal_total = rowSums(combined[,c("Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970") ], na.rm = TRUE), 
         Municipal_eligible = rowSums(!is.na(combined[,c("Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970") ])), 
         EU_total = rowSums(combined[,c("EU_2019", "EU_2009")], na.rm = TRUE), 
         EU_eligible = rowSums(!is.na(combined[,c("EU_2019", "EU_2009") ])), 
         Total_total = rowSums(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970", 
                                                   "Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970", 
                                                   "Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970", 
                                                   "EU_2019", "EU_2009")], na.rm = TRUE), 
         Total_eligible = rowSums(!is.na(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970", 
                                                                "Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970", 
                                                                "Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970", 
                                                                "EU_2019", "EU_2009")])), 
         Domestic_total = rowSums(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970", 
                                                   "Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970", 
                                                   "Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970")], na.rm = TRUE), 
         Domestic_eligible = rowSums(!is.na(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970", 
                                                                "Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970", 
                                                                "Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970")])),
                                       National_prop = National_total/National_eligible, 
                                       Regional_prop = Regional_total/Regional_eligible, 
                                       Municipal_prop = Municipal_total/Municipal_eligible, 
                                       EU_prop = EU_total/EU_eligible, 
         Total_prop = Total_total/Total_eligible, 
         Domestic_prop = Domestic_total/Domestic_eligible, 
         Domestic_scaled = scale(Domestic_prop), 
         Total_scaled = scale(Total_prop), 
         Regional_scaled = scale(Regional_prop), 
         Municipal_scaled = scale(Municipal_prop), 
         National_scaled = scale(National_prop), 
         EU_scaled = scale(EU_prop))
                          
#Filtering for complete sets

pairs = combined %>% 
  group_by(LopNrParID) %>% 
  filter(is.na(trust) == F & is.na(behavioral) ==F) %>% 
  filter( n() > 1 ) %>%
  mutate(TWINID = rep(c(1, 2),length.out = n())) %>% 
  ungroup()

#Stacking data by twin pairs

pairs = pairs %>% 
  select(LopNr, LopNrParID, trust, oldtrust, trust_scaled, behavioral, contactpol, contactpub, protest, boycott, finance, petition, TWINID, BESTZYG, SEX, Byear, National_prop, Regional_prop, Municipal_prop, Domestic_prop, Total_prop, interpersonal_1, interpersonal_2, interpersonal_scaled) 

pairs$current_date <- as.numeric(format(Sys.Date(), "%Y%m"))
pairs$current_year <- floor(pairs$current_date/100)
pairs$current_month <- pairs$current_date %% 100

# Calculate age in years and adjust for month of birth
pairs$Byear = as.numeric(pairs$Byear)
pairs$age <- pairs$current_year - floor(pairs$Byear/100) - 1
pairs$age = ifelse(pairs$current_month >= (pairs$Byear %% 100), pairs$age <- pairs$age + 1, pairs$age <- pairs$age)

twin1 = pairs %>% 
  filter(TWINID == 1)

twin2 = pairs %>% 
  filter(TWINID == 2)

final = merge(twin1, twin2, by.x = "LopNrParID", by.y = "LopNrParID", suffixes = c("_T1", "_T2"))

#Creating MZ and DZ datasets
SALTY_MZ <- subset(final, BESTZYG_T1 == 1)
SALTY_DZ <- subset(final, BESTZYG_T1 != 1)

```

```{r Bivariate Cholesky behavioral scale} 

CHO1 <- umxACE(selDVs = c("trust_scaled", "behavioral"), selCovs = c("age", "SEX"), sep = "_T", mzData = SALTY_MZ, dzData = SALTY_DZ, tryHard = "yes", optimizer = "SLSQP")

genCorAlg <- mxAlgebra(cov2cor(top.a%*%t(top.a)), name="geneticCorrelation")
genCorCI <- mxCI("geneticCorrelation")

comCorAlg <- mxAlgebra(cov2cor(top.c%*%t(top.c)), name="comCorrelation")
comCorCI <- mxCI("comCorrelation")

envCorAlg <- mxAlgebra(cov2cor(top.e%*%t(top.e)), name="envCorrelation")
envCorCI <- mxCI("envCorrelation")

CHO1 <- mxModel(CHO1, genCorAlg, genCorCI, envCorAlg, envCorCI, comCorAlg, comCorCI)

CHO1 <- mxRun(CHO1, intervals = T)

summary(CHO1, verbose = T)
```

```{r Proportional genetic correlation}

genCor = CHO1$geneticCorrelation$result[1,2]

gen_var_behav = CHO1$output$algebras$top.A_std[1,1]
gen_var_trust = CHO1$output$algebras$top.A_std[2,2]

gen_cov <- genCor * sqrt(gen_var_behav) * sqrt(gen_var_trust)


```

```{r Proportional common correlation}

comCor = CHO1$comCorrelation$result[2,2]

com_var_behav = CHO1$output$algebras$top.C_std[1,1]
com_var_trust = CHO1$output$algebras$top.C_std[2,2]

com_cov <- comCor * sqrt(com_var_behav) * sqrt(com_var_trust)

```

```{r Proportional environmental correlation}

envCor = CHO1$envCorrelation$result[1,2]

env_var_behav = CHO1$output$algebras$top.E_std[1,1]
env_var_trust = CHO1$output$algebras$top.E_std[2,2]

env_cov <- envCor * sqrt(env_var_behav) * sqrt(env_var_trust)

gen_perc <- (gen_cov / (env_cov + gen_cov + com_cov)) * 100

com_perc <- (com_cov / (env_cov + gen_cov + com_cov)) * 100

env_perc <- (env_cov / (env_cov + gen_cov + com_cov)) * 100

gen_perc

env_perc

com_perc

```

```{r Confidence intervals, genetic correlation, behavioral}
require(boot)

my_boot <- function(data, indices) {
  SALTY_MZ <- subset(data[indices,], BESTZYG_T1 == 1)
  SALTY_DZ <- subset(data[indices,], BESTZYG_T1 != 1)
  
  # run the Cholesky model with umxACE
  CHO1 <- umxACE(selDVs = c("trust_scaled", "behavioral"), selCovs = c("age", "SEX"), sep = "_T", mzData = SALTY_MZ, dzData = SALTY_DZ, tryHard = "yes", optimizer = "SLSQP")
  
  CHO1 = umxConfint(CHO1, optimizer = "SLSQP", run = T)

  genCorAlg <- mxAlgebra(cov2cor(top.a%*%t(top.a)), name="geneticCorrelation")
  genCorCI <- mxCI("geneticCorrelation")
  
  comCorAlg <- mxAlgebra(cov2cor(top.c%*%t(top.c)), name="comCorrelation")
  comCorCI <- mxCI("comCorrelation")

  envCorAlg <- mxAlgebra(cov2cor(top.e%*%t(top.e)), name="envCorrelation")
  envCorCI <- mxCI("envCorrelation")

  CHO1 <- mxModel(CHO1, genCorAlg, genCorCI, envCorAlg, envCorCI, comCorAlg, comCorCI)

  CHO1 <- mxRun(CHO1, intervals = T)
  
  # compute the genetic covariance
  genCor <- CHO1$geneticCorrelation$result[1,2]
  gen_var_behav <- CHO1$output$algebras$top.A_std[1,1]
  gen_var_trust <- CHO1$output$algebras$top.A_std[2,2]
  gen_cov <- genCor * sqrt(gen_var_behav) * sqrt(gen_var_trust)
  
  comCor = CHO1$comCorrelation$result[2,2]
  com_var_behav = CHO1$output$algebras$top.C_std[1,1]
  com_var_trust = CHO1$output$algebras$top.C_std[2,2]
  com_cov <- comCor * sqrt(com_var_behav) * sqrt(com_var_trust)
  
  
  envCor = CHO1$envCorrelation$result[1,2]
  env_var_behav = CHO1$output$algebras$top.E_std[1,1]
  env_var_trust = CHO1$output$algebras$top.E_std[2,2]
  env_cov <- envCor * sqrt(env_var_behav) * sqrt(env_var_trust)
  # compute the proportion of the phenotypic correlation attributable to genetic factors
  gen_perc <- (gen_cov / (env_cov + gen_cov + com_cov)) * 100
  
  return(gen_perc)
}

# run boot() with my_boot as the statistic function
set.seed(123)
genprop_behav <- boot(data = final, statistic = my_boot, R = 5)
genprop_behav_ci <- boot.ci(genprop_behav, type = "basic", conf = 0.95)
genprop_behav
genprop_behav_ci
```

```{r Confidence intervals, common correlation, behavioral}
require(boot)

my_boot <- function(data, indices) {
  SALTY_MZ <- subset(data[indices,], BESTZYG_T1 == 1)
  SALTY_DZ <- subset(data[indices,], BESTZYG_T1 != 1)
  
  # run the Cholesky model with umxACE
  CHO1 <- umxACE(selDVs = c("trust", "behavioral"), selCovs = c("age", "SEX"), sep = "_T", mzData = SALTY_MZ, dzData = SALTY_DZ, tryHard = "yes", optimizer = "SLSQP")
  
  CHO1 = umxConfint(CHO1, optimizer = "SLSQP", run = T)
  
  genCorAlg <- mxAlgebra(cov2cor(top.a%*%t(top.a)), name="geneticCorrelation")
  genCorCI <- mxCI("geneticCorrelation")
  
  comCorAlg <- mxAlgebra(cov2cor(top.c%*%t(top.c)), name="comCorrelation")
  comCorCI <- mxCI("comCorrelation")

  envCorAlg <- mxAlgebra(cov2cor(top.e%*%t(top.e)), name="envCorrelation")
  envCorCI <- mxCI("envCorrelation")

  CHO1 <- mxModel(CHO1, genCorAlg, genCorCI, envCorAlg, envCorCI, comCorAlg, comCorCI)

  CHO1 <- mxRun(CHO1, intervals = T)
  
  # compute the genetic covariance
  genCor <- CHO1$geneticCorrelation$result[1,2]
  gen_var_behav <- CHO1$output$algebras$top.A_std[1,1]
  gen_var_trust <- CHO1$output$algebras$top.A_std[2,2]
  gen_cov <- genCor * sqrt(gen_var_behav) * sqrt(gen_var_trust)
  
  comCor = CHO1$comCorrelation$result[2,2]
  com_var_behav = CHO1$output$algebras$top.C_std[1,1]
  com_var_trust = CHO1$output$algebras$top.C_std[2,2]
  com_cov <- comCor * sqrt(com_var_behav) * sqrt(com_var_trust)
  
  
  envCor = CHO1$envCorrelation$result[1,2]
  env_var_behav = CHO1$output$algebras$top.E_std[1,1]
  env_var_trust = CHO1$output$algebras$top.E_std[2,2]
  env_cov <- envCor * sqrt(env_var_behav) * sqrt(env_var_trust)
  
  # compute the proportion of the phenotypic correlation attributable to genetic factors
  com_perc <- (com_cov / (env_cov + gen_cov + com_cov)) * 100
  
  return(com_perc)
}

# run boot() with my_boot as the statistic function
set.seed(123)
comprop_behav <- boot(data = final, statistic = my_boot, R = 5)
comprop_behav_ci <- boot.ci(comprop_behav, type = "norm", conf = 0.95)
comprop_behav
comprop_behav_ci
```

```{r Confidence intervals, environmental correlation, behavioral}
require(boot)

my_boot <- function(data, indices) {
  SALTY_MZ <- subset(data[indices,], BESTZYG_T1 == 1)
  SALTY_DZ <- subset(data[indices,], BESTZYG_T1 != 1)
  
  # run the Cholesky model with umxACE
  CHO1 <- umxACE(selDVs = c("trust_scaled", "behavioral"), selCovs = c("age", "SEX"), sep = "_T", mzData = SALTY_MZ, dzData = SALTY_DZ, tryHard = "yes", optimizer = "SLSQP")
  
  CHO1 = umxConfint(CHO1, optimizer = "SLSQP", run = T)

  genCorAlg <- mxAlgebra(cov2cor(top.a%*%t(top.a)), name="geneticCorrelation")
  genCorCI <- mxCI("geneticCorrelation")
  
  comCorAlg <- mxAlgebra(cov2cor(top.c%*%t(top.c)), name="comCorrelation")
  comCorCI <- mxCI("comCorrelation")

  envCorAlg <- mxAlgebra(cov2cor(top.e%*%t(top.e)), name="envCorrelation")
  envCorCI <- mxCI("envCorrelation")

  CHO1 <- mxModel(CHO1, genCorAlg, genCorCI, envCorAlg, envCorCI, comCorAlg, comCorCI)

  CHO1 <- mxRun(CHO1, intervals = T)
  
  # compute the genetic covariance
  genCor <- CHO1$geneticCorrelation$result[1,2]
  gen_var_behav <- CHO1$output$algebras$top.A_std[1,1]
  gen_var_trust <- CHO1$output$algebras$top.A_std[2,2]
  gen_cov <- genCor * sqrt(gen_var_behav) * sqrt(gen_var_trust)
  
  comCor = CHO1$comCorrelation$result[2,2]
  com_var_behav = CHO1$output$algebras$top.C_std[1,1]
  com_var_trust = CHO1$output$algebras$top.C_std[2,2]
  com_cov <- comCor * sqrt(com_var_behav) * sqrt(com_var_trust)
  
  
  envCor = CHO1$envCorrelation$result[1,2]
  env_var_behav = CHO1$output$algebras$top.E_std[1,1]
  env_var_trust = CHO1$output$algebras$top.E_std[2,2]
  env_cov <- envCor * sqrt(env_var_behav) * sqrt(env_var_trust)
  
  # compute the proportion of the phenotypic correlation attributable to genetic factors
  env_perc <- (env_cov / (env_cov + com_cov + gen_cov)) * 100
  
  return(env_perc)
}

# run boot() with my_boot as the statistic function
set.seed(123)
envprop_behav <- boot(data = final, statistic = my_boot, R = 5)
envprop_behav_ci <- boot.ci(envprop_behav, type = "basic", conf = 0.95)
envprop_behav
envprop_behav_ci
```

```{r Bivariate Cholesky turnout scale} 

CHO1 <- umxACE(selDVs = c("trust_scaled", "National_prop"), selCovs = c("age", "SEX"), sep = "_T", mzData = SALTY_MZ, dzData = SALTY_DZ, tryHard = "no", optimizer = "SLSQP")

genCorAlg <- mxAlgebra(cov2cor(top.a%*%t(top.a)), name="geneticCorrelation")
genCorCI <- mxCI("geneticCorrelation")

comCorAlg <- mxAlgebra(cov2cor(top.c%*%t(top.c)), name="comCorrelation")
comCorCI <- mxCI("comCorrelation")

envCorAlg <- mxAlgebra(cov2cor(top.e%*%t(top.e)), name="envCorrelation")
envCorCI <- mxCI("envCorrelation")

CHO1 <- mxModel(CHO1, genCorAlg, genCorCI, envCorAlg, envCorCI, comCorAlg, comCorCI)

CHO1 <- mxRun(CHO1, intervals = T)

summary(CHO1, verbose = T)
```
